##
## The downloaded binary packages are in
## /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpQoXYeb/downloaded_packages
##
## The downloaded binary packages are in
## /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpQoXYeb/downloaded_packages
##
## The downloaded binary packages are in
## /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpQoXYeb/downloaded_packages
SO1 <- SCTransform(SO) %>%
RunPCA() %>%
FindNeighbors(dims = 1:30) %>%
FindClusters() %>%
RunUMAP(dims = 1:30)## Running SCTransform on assay: RNA
## Warning: The `slot` argument of `GetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
## Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
## Calculating cell attributes from input UMI matrix: log_umi
## Variance stabilizing transformation of count matrix of size 21397 by 5010
## Model formula is y ~ log_umi
## Get Negative Binomial regression parameters per gene
## Using 2000 genes, 5000 cells
## Found 126 outliers - those will be ignored in fitting/regularization step
## Second step: Get residuals using fitted parameters for 21397 genes
## Computing corrected count matrix for 21397 genes
## Calculating gene attributes
## Wall clock passed: Time difference of 11.76 secs
## Determine variable features
## Centering data matrix
## Place corrected count matrix in counts slot
## Warning: The `slot` argument of `SetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
## Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Different cells and/or features from existing assay SCT
## Set default assay to SCT
## PC_ 1
## Positive: Pappa2, Malat1, S100g, Zfand5, Pde10a, Ramp3, Robo2, Sdc4, Itga4, Gm37376
## Aard, Nadk2, Nos1, Ranbp3l, Wwc2, Sgms2, Neat1, Irx1, Col4a3, Col4a4
## Tmem158, Aebp1, Cdkn1c, Itprid2, Cacna1d, Ptgs2, Ccdc80, Bmp3, Camk2d, Uroc1
## Negative: Umod, Egf, Tmem52b, Fabp3, Klk1, Sostdc1, Sult1d1, Prdx5, Wfdc15b, Ly6a
## Foxq1, Cox6c, Krt7, Ckb, Cldn19, Cox5b, Slc25a5, Wfdc2, Atp5g1, Atp1a1
## Ggt1, Cox4i1, Atp1b1, Gm47708, Gadd45g, Cox7b, Ndufa4, Cox8a, Atp5b, Gm44120
## PC_ 2
## Positive: Ubb, Mt1, Fth1, Gm22133, Hspa1a, Ldhb, Actb, Hspa1b, H3f3b, Fos
## Junb, Ftl1, Eif1, Jund, Jun, Cd63, Fxyd2, Prdx1, Mt2, Tmem213
## Mpc2, Cd9, Mgst1, Rps24, Rpl37, Ier2, Cox8a, Btg2, Gm8355, Egr1
## Negative: Mir6236, Gm37376, CT010467.1, Egf, Malat1, Umod, Gm24447, Nme7, mt-Nd5, mt-Rnr1
## Slc12a1, mt-Nd6, Kcnq1ot1, Neat1, Tmem52b, Lars2, Atp1b1, mt-Nd4l, Dst, Etnk1
## Wnk1, mt-Rnr2, Atrx, Slc5a3, Syne2, Robo2, mt-Co1, mt-Nd2, Zbtb20, Pnisr
## PC_ 3
## Positive: Gm22133, Ldhb, Mgst1, Car15, Mpc2, Cd63, Tmem213, Fxyd2, Cd9, Cldn10
## Ftl1, Rpl9, Prdx1, Rplp1, Wfdc2, Rpl7, S100a1, Fth1, Ramp3, Rps27a
## Rps24, Ly6a, Irx1, Ppp1r1a, Atp5e, Apoe, Mdh1, Cox8a, Tmbim6, Rpl37
## Negative: Fos, Jun, Junb, Egr1, Hspa1b, Btg2, Hspa1a, Atf3, Zfp36, Klf6
## Ier2, Fosb, Socs3, Klf2, Jund, Dnajb1, Gadd45g, Tsc22d1, Gm37376, Rhob
## Malat1, Gadd45b, Actb, H1f2, Ubc, H2bc4, H3f3b, Egf, Nr4a1, Dusp1
## PC_ 4
## Positive: S100g, Gm8420, Actb, Sdc4, Pth1r, Igfbp7, Gm8885, Pappa2, Malat1, Umod
## Ramp3, Tmem52b, Egf, Cebpd, Cd9, Tmsb10, Sgms2, Wnk1, Gm28037, Rpl41
## Ppm1h, Tmem158, Ccdc80, Lmna, Tbxas1, Pdcd4, Ppp1r1a, Tmbim6, Col4a3, Gm11808
## Negative: Mt1, Mt2, Rpl15-ps2, Rpl13, Ptger3, CT010467.1, Rpl10, Rpl9, Gm22133, mt-Rnr2
## Apoe, Rpl24, Aebp1, Rpl17, Jun, Rpsa, Ckb, Hspa8, Eef1a1, Rpl18
## Cd63, Gm13340, Gpx6, Gm5905, Uba52-ps, Gm23935, Gpx4, Rplp1, Mfsd4a, Chchd10
## PC_ 5
## Positive: Pappa2, Tmem52b, Umod, Aard, Egf, Sult1d1, Foxq1, Actb, Tmem158, Gm44120
## Mcub, Pde10a, Rpl15-ps2, Defb42, Tmsb4x, Dctd, Car15, Rpl10, Kcnj10, 5330417C22Rik
## Begain, Rpl13, Rpl24, Abca13, Gsn, Zfand5, Eef1a1, Fabp3, 1700028P14Rik, Gm22133
## Negative: mt-Nd5, mt-Co1, mt-Nd4l, Gm28437, mt-Cytb, Malat1, Gm28438, mt-Rnr1, Kng1, Nme7
## mt-Nd2, Mir6236, Gm8420, Hspa1b, Gm28037, Gm10925, Defb1, Rps21, Atp1b1, Fxyd2
## Klk1, Gm11808, Gm8885, Tnfaip2, Cox6a1, Gm12338, S100g, Sostdc1, mt-Nd3, Car12
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 5010
## Number of edges: 184583
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7657
## Number of communities: 10
## Elapsed time: 0 seconds
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 11:42:13 UMAP embedding parameters a = 0.9922 b = 1.112
## 11:42:13 Read 5010 rows and found 30 numeric columns
## 11:42:13 Using Annoy for neighbor search, n_neighbors = 30
## 11:42:13 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 11:42:13 Writing NN index file to temp file /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpQoXYeb/file157107176d16a
## 11:42:13 Searching Annoy index using 1 thread, search_k = 3000
## 11:42:14 Annoy recall = 100%
## 11:42:14 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 11:42:14 Initializing from normalized Laplacian + noise (using RSpectra)
## 11:42:14 Commencing optimization for 500 epochs, with 212904 positive edges
## 11:42:14 Using rng type: pcg
## 11:42:17 Optimization finished
## Calculating cluster SO1
## Warning: `PackageCheck()` was deprecated in SeuratObject 5.0.0.
## ℹ Please use `rlang::check_installed()` instead.
## ℹ The deprecated feature was likely used in the Seurat package.
## Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
## (default method for FindMarkers) please install the presto package
## --------------------------------------------
## install.packages('devtools')
## devtools::install_github('immunogenomics/presto')
## --------------------------------------------
## After installation of presto, Seurat will automatically use the more
## efficient implementation (no further action necessary).
## This message will be shown once per session
## Calculating cluster SO4
SO1.markers %>%
group_by(cluster) %>%
dplyr::filter(avg_log2FC > 1) %>%
slice_head(n = 10) %>%
ungroup() -> top10
DoHeatmap(SO1, features = top10$gene) + NoLegend()## Warning in DoHeatmap(SO1, features = top10$gene): The following features were
## omitted as they were not found in the scale.data slot for the SCT assay:
## Ldhb-ps, Gm8292, Rps6
## Warning: The `slot` argument of `FetchData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
## Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
markers.to.plot1 <- c("Ptger3","Hspa8","S100g")
DotPlot(SO1,
features = markers.to.plot1,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()SO1 expresses highly in Ptger3, SO4 expresses high in S100g.
S01 Also has a lot more Ribosomal gene influence.
Why?
df<- SO1_markers %>% arrange(desc(avg_log2FC))
df2 <- df %>% filter(p_val_adj < 0.05)
DEG_list <- df2
markers <- DEG_list %>% rownames_to_column(var="SYMBOL")
head(markers, n = 50)ENTREZ_list <- bitr(
geneID = rownames(DEG_list),
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Mm.eg.db
)## 'select()' returned 1:1 mapping between keys and columns
## Warning in bitr(geneID = rownames(DEG_list), fromType = "SYMBOL", toType =
## "ENTREZID", : 4.37% of input gene IDs are fail to map...
markers <- ENTREZ_list %>% inner_join(markers, by = "SYMBOL")
# Removing genes that are not statistically significant.
markers <- markers %>% dplyr::filter(p_val_adj < 0.05)
#head(markers, n = 50)
pos.markers <- markers %>% dplyr::filter(avg_log2FC > 0) %>% arrange(desc(abs(avg_log2FC)))
#head(pos.markers, n = 50)
pos.ranks <- pos.markers$ENTREZID[abs(pos.markers$avg_log2FC) > 0]
#head(pos.ranks)
pos_go <- enrichGO(gene = pos.ranks, #a vector of entrez gene id
OrgDb = "org.Mm.eg.db",
ont = "BP",
readable = TRUE) #whether mapping gene ID to gene Name
pos_go## #
## # over-representation test
## #
## #...@organism Mus musculus
## #...@ontology BP
## #...@keytype ENTREZID
## #...@gene chr [1:1619] "20564" "380785" "94283" "15080" "54377" "622570" "69524" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...1028 enriched terms found
## 'data.frame': 1028 obs. of 12 variables:
## $ ID : chr "GO:0002181" "GO:0022613" "GO:0006457" "GO:0008380" ...
## $ Description : chr "cytoplasmic translation" "ribonucleoprotein complex biogenesis" "protein folding" "RNA splicing" ...
## $ GeneRatio : chr "58/1431" "92/1431" "49/1431" "79/1431" ...
## $ BgRatio : chr "171/28928" "452/28928" "185/28928" "464/28928" ...
## $ RichFactor : num 0.339 0.204 0.265 0.17 0.261 ...
## $ FoldEnrichment: num 6.86 4.11 5.35 3.44 5.27 ...
## $ zScore : num 17.5 15.2 13.6 12.1 13.3 ...
## $ pvalue : num 6.44e-33 1.27e-31 1.28e-22 4.21e-22 7.06e-22 ...
## $ p.adjust : num 3.66e-29 3.63e-28 2.42e-19 5.99e-19 8.03e-19 ...
## $ qvalue : num 2.65e-29 2.62e-28 1.75e-19 4.33e-19 5.81e-19 ...
## $ geneID : chr "Gm6133/Rps6/Rpl13/Rpl18/Eif2s3y/Rpl17/Rpsa/Rps19/Rpl9/Rpl10a/Rpl24/Rpl28/Rpl18a/Rpl21/Rpl27a/Rps6-ps4/Rpl12/Rpl"| __truncated__ "Rps6/Eif2s3y/Nol3/Mrto4/Bysl/Rpsa/Nop56/Celf2/Rps19/Rpl24/Utp14b/Fbl/Zfp593/Gemin7/Mphosph10/Rps6-ps4/Srsf6/Fts"| __truncated__ "Hspa8/Fkbp11/Ric3/Umod/Hspb1/Hsph1/Ranbp2/Dnajc2/Dnajc10/Ahsa2/Tor2a/Nktr/Ppil1/Mkks/Hspa4/Chordc1/Cct3/Aip/Dna"| __truncated__ "Hspa8/Nol3/Celf2/Tsen34/U2af1/Umod/Srsf3/Polr2a/Gemin7/Srsf6/Srsf7/Zc3h13/Ybx1/Mfap1a/Prpf8/Zfp326/Rnpc3/Phf5a/"| __truncated__ ...
## $ Count : int 58 92 49 79 48 63 48 74 64 23 ...
## #...Citation
## G Yu. Thirteen years of clusterProfiler. The Innovation. 2024, 5(6):100722
dotplot(pos_go) +
ggtitle("") +
theme_classic() +
theme(
plot.title = element_text(hjust = 0.5),
legend.position = "left",
axis.text.y = element_text(hjust = 0, size = 10)) +
scale_y_discrete(position = "right",
labels = function(x) str_wrap(x, width = 25)) # Wrap y-axis labels to 2 lines## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
I Think that SO1, has too much ribosomal influence so I will try to remove it here and see what the changes are.
df<- SO4_markers %>% arrange(desc(avg_log2FC))
df2 <- df %>% filter(p_val_adj < 0.05)
DEG_list <- df2
markers <- DEG_list %>% rownames_to_column(var="SYMBOL")
head(markers, n = 50)ENTREZ_list <- bitr(
geneID = rownames(DEG_list),
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Mm.eg.db
)## 'select()' returned 1:1 mapping between keys and columns
## Warning in bitr(geneID = rownames(DEG_list), fromType = "SYMBOL", toType =
## "ENTREZID", : 8.82% of input gene IDs are fail to map...
markers <- ENTREZ_list %>% inner_join(markers, by = "SYMBOL")
# Removing genes that are not statistically significant.
markers <- markers %>% dplyr::filter(p_val_adj < 0.05)
#head(markers, n = 50)
pos.markers <- markers %>% dplyr::filter(avg_log2FC > 0) %>% arrange(desc(abs(avg_log2FC)))
#head(pos.markers, n = 50)
pos.ranks <- pos.markers$ENTREZID[abs(pos.markers$avg_log2FC) > 0]
#head(pos.ranks)
pos_go <- enrichGO(gene = pos.ranks, #a vector of entrez gene id
OrgDb = "org.Mm.eg.db",
ont = "BP",
readable = TRUE) #whether mapping gene ID to gene Name
pos_go## #
## # over-representation test
## #
## #...@organism Mus musculus
## #...@ontology BP
## #...@keytype ENTREZID
## #...@gene chr [1:424] "667014" "100041678" "667937" "671215" "13653" "384179" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...377 enriched terms found
## 'data.frame': 377 obs. of 12 variables:
## $ ID : chr "GO:0140236" "GO:0140241" "GO:0140242" "GO:0002181" ...
## $ Description : chr "translation at presynapse" "translation at synapse" "translation at postsynapse" "cytoplasmic translation" ...
## $ GeneRatio : chr "14/306" "14/306" "14/306" "21/306" ...
## $ BgRatio : chr "47/28928" "48/28928" "48/28928" "171/28928" ...
## $ RichFactor : num 0.2979 0.2917 0.2917 0.1228 0.0806 ...
## $ FoldEnrichment: num 28.16 27.57 27.57 11.61 7.62 ...
## $ zScore : num 19.27 19.05 19.05 14.39 9.37 ...
## $ pvalue : num 4.07e-17 5.69e-17 5.69e-17 1.67e-16 1.54e-09 ...
## $ p.adjust : num 7.33e-14 7.33e-14 7.33e-14 1.61e-13 1.19e-06 ...
## $ qvalue : num 5.46e-14 5.46e-14 5.46e-14 1.20e-13 8.85e-07 ...
## $ geneID : chr "Rpl38/Rps28/Rpl37a/Rpl37/Rps26/Rps27/Rpl23a/Rpl36a/Rpl34/Rps27a/Rpl36/Rps12/Rpl22/Rpl35a" "Rpl38/Rps28/Rpl37a/Rpl37/Rps26/Rps27/Rpl23a/Rpl36a/Rpl34/Rps27a/Rpl36/Rps12/Rpl22/Rpl35a" "Rpl38/Rps28/Rpl37a/Rpl37/Rps26/Rps27/Rpl23a/Rpl36a/Rpl34/Rps27a/Rpl36/Rps12/Rpl22/Rpl35a" "Rpl38/Rps28/Rpl37a/Rpl37/Rpl41/Rpl39/Rps29/Rps21/Rpl31/Pabpc1/Rps26/Rpl22l1/Rpl6l/Rpl23a/Rpl36a/Rpl34/Rps27a/Rp"| __truncated__ ...
## $ Count : int 14 14 14 21 15 16 21 19 18 19 ...
## #...Citation
## G Yu. Thirteen years of clusterProfiler. The Innovation. 2024, 5(6):100722
dotplot(pos_go) +
ggtitle("") +
theme_classic() +
theme(
plot.title = element_text(hjust = 0.5),
legend.position = "left",
axis.text.y = element_text(hjust = 0, size = 10)) +
scale_y_discrete(position = "right",
labels = function(x) str_wrap(x, width = 25)) # Wrap y-axis labels to 2 lines## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.